An Efficient Geometric Integrator for Thermostatted 
Anti-/Ferromagnetic Models 



Teijo Arponen Ben Leimkuhler 

Institute of Mathematics Department of Mathematics 
Helsinki University of Technology University of Leicester 

Finland U.K. 

1st February 2008 



Abstract 

(Anti)-/ferromagnetic Heisenberg spin models arise from discretization of Landau- 
Lifshitz models in micromagnetic modelling. In many applications it is essential to 
study the behavior of the system at a fixed temperature. A formulation for ther- 
mostatted spin dynamics was given by Bulgac and Kusnetsov [Sj , which incorporates 
a complicated nonlinear dissipation/driving term while preserving spin length. It 
is essential to properly model this term in simulation, and simplified schemes give 
poor numerical performance, e.g. requiring an excessively small timestep for stable 
integration. In this paper we present an efficient, structure-preserving method for 
thermostatted spin dynamics. 

Keywords: Heisenberg ferromagnet, micromagnetics, spin dynamics, Landau-Lifschitz 
equation, Gilbert damping, thermostats, constant temperature, domain walls, geometric 
integrator, reversible method 

1 Introduction 

In recent years geometric integrators have become ubiquitous for numerical treatment of 
differential equations. By a geometric integrator is meant a numerical method that pre- 
serves some known structure of the continuous flow. Geometric integrators are particularly 
important for long term simulations, as used in molecular sampling or celestial mechanics. 
In this paper we consider the application of geometric integration principles for the types 
of spin dynamics systems which arise frequently in modelling of ferromagnets and anti- 
ferromagnets. Efficient Lie-Poisson schemes for classical spin dynamics described by the 
Landau-Lifshitz (LL) equation were studied in ^U], and related multisymplectic schemes in 
10]. Here we develop and test a geometric integrator for a semi-discrete Landau-Lifshitz- 
Gilbert (LLG) equation which includes a nonlinear dissipative term. This dissipative 
system forms the foundation for a more complicated thermostatted model, following the 
approach of Bulgac and Kusnetsov I3I2]- We design an effective splitting technique for 
the full coupled system. 

LL and LLG are currently a very active topic of research. Other approaches to them 
can be found in [HI ESI El HDl UHl 1211 HH HZ] , who also provide further references. However 
none of these consider a thermostatted version. 
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Simulation with the thermostatted version shows fascinating global behavior: the 
system first arranges into patterns (spin domains) with slowly moving domain walls, then 
goes into a quasi-chaotic state and quickly rearranges itself into completely new spin 
domains. This kind of transition would not be possible with local interactions only. Here 
the thermostatting variable is defined in such a way that it has a global character. 

A number of recent articles have focussed on the geometric integration of molecular 
systems in the canonical ensemble |31 El El E]- In these articles, the aim has been to 
start from a Hamiltonian formulation for thermostatted molecular simulation and then 
to provide a suitable symplectic integrator. The starting point is usually Nose dynamics, 
although generalizations are possible. 

Since the constant energy Heisenberg spin system is Lie-Poisson, it is natural to seek 
a Lie-Poisson system to model the action of the thermostat. While it is possible (with 
some additional complication, due to the presence of constraints) to develop such a model 
for the thermostatted Heisenberg model, based on the ideas in ^2j, it is much different in 
character from the corresponding molecular dynamics models (see the appendix). In par- 
ticular, this approach appears to require introduction of many thermostatting variables 
which act differently on each spin vector of the system. In the context of magnetic mod- 
els, this approach therefore sacrifices an important feature of Nose molecular dynamics: 
the apparent compatibility between the thermostatted quasi-dynamics and the micro- 
canonical dynamics. (Even though Nose dynamics is typically only validated based on 
a phase- sampling correspondence, there is widespread agreement that the thermostatted 
dynamics is relevant for modelling dynamics of an appropriate extended system which is 
not too different in character from the microcanonical version.) Moreover, the Lie-Poisson 
thermostats add additional complexity in the form of a relatively complex bath model. 

Given these complications, we believe the best available starting point for geometric 
integration of thermostatted spin dynamics is the alternative framework of Bulgac and 
Kusnetsov, based loosely on Nose-Hoover (NH) dynamics. Like NH molecular dynamics, 
these formulations sacrifice Hamiltonian structure, while retaining a reversing symmetry. 
It is unclear the extent to which this loss of structure affects the stability of methods and 
the ultimate resolution of macroscopic features of the spin model. Although in molecular 
dynamics it is known that the reversible-only methods are often inferior to their symplectic 
counterparts [13], it is also well established that NH-type methods are far superior to 
methods that are neither symplectic nor reversible. 

The rest of the paper is organized as follows: in Section 101 we review splitting methods 
and apply them to our models. In Sections 1210] we present the models and methods in 
detail. In Section El we present numerical results. Finally in Section IHl we present some 
conclusions and discussion. 

1.1 Background: Review of splitting methods 

The reader is referred to fH^ ITH] for a detailed discussion of splitting methods. To briefly 
describe their basic construction, consider a differential equation u = f{u), with flow map 
If / = /i + /2, we have = ^rji <=> ^rj2 + 0{t'^). If the flows on vector fields /, 
/i, and /2 share a first integral, then the composed map will preserve it as well. In this 
way, geometric integrators can be developed to preserve general classes of Lie groups. If 
the vector field is time-reversible, i.e. f{Ru) = —Rf{u) for some linear involution then 
the symmetric concatenation or "Strang Splitting" = ^^T,h ° ^tJ2 ° ^^T,h' where 
/i, /2 are reversible vector fields, gives a time-reversible map (i?<|-i = d^joi?), which, 
moreover, provides a second-order approximation of the solution on a finite time interval. 



As an example, ii H = H{q,p) = T{j>) + V{q), the leapfrog (Stormer/Verlet) integrator 
results from the concatenation ^^-u = $1^ v o^tT ° ^i^v- 

' 2, ' ' 2 ' 

The construction of splitting methods for various types of flows, and with various orders 
of accuracy, is discussed in a number of papers (see, e.g, [211 ED])- Practical splitting-based 
geometric integrators have been constructed by mathematicians, chemists and physicists 
for a wide variety of important applications, including the rigid body, general holonomic 
constraints, particle accelerator models, and the solar system. Vector field splittings were 
used in ^UI to obtain efficient time-reversible integrators for (undamped) spin systems; it 
is this fundamental scheme that we have extended in this paper to treatment of dissipative 
and thermostatted systems. 



2 The original Landau-Lifshitz model as a Poisson 
system 

There are several versions of the Landau-Lifshitz equation depending on which forces and 
fields are taken into account. The version we use here is that of |7|j, discarding the external 
and demagnetizing field. (Schemes for more general formulations would build on the work 
presented here.) The equation can be written in the form: 

^^S = S yiV^S + S DS, (1) 

where xG/x/c]R^,/an interval, S{x,t) is a unit vector in representing the classical 
spin at position x and time t, and D is a diagonal matrix representing anisotropy. Clearly 
IS"! = constant in time: 

-|^(x,t)|2 = 2S{x,t) ■ g-^S{x,t) = Vx,t. 

Following the usual practice, we discretize the spatial variable x using second order central 
differences on a regular lattice as in so that in the discretized system the unit length 
property is conserved. We then get a Poisson system on a lattice. Without loss of 
generality we may assume the lattice size to be 1: 

S (x, ■) I > Zij 
V^5'(x, ■) (-^ Zij-i + Zij+i + Zi-ij + Zi+ij — Azij, 

hence (P) becomes 

Zij = Zij X (^Zij — i -\- Zij^i -\- Zi—ij -\- Zi^ij ^Zij) -\- Zij X DZij. (2) 

Note that the —^Zij term can be dropped out. Here we have an n x n lattice of spins: the 
variable Zij is on the unit sphere of when i,j G {1, . . . When either i or j index 
is zero or n -|- 1, those represent boundaries. Except for the case of periodic boundary 
conditions, these boundary terms are different from the spins: they are an artefact of 
discretization, and do not have a counterpart in the continuum case (^. Especially, they 
are not necessarily of unit length. We do not represent equations of motion to them, hence 
they are assumed constants. 

By periodic boundary conditions we mean 



Zoj Ziij^ ZiQ Zin- 



(3) 



Next we define the Poisson structure matrix. Let us denote 

._ \ T \ T \ \ T \ T \ T \ I l'^ 

^ L^ll I ^12 I • • • I ^In I ^21 I ^22 I • • • I ^nn\ ' 

i.e. 2; is a column vector. For an arbitrary v =: [a, 6, c]"^ G we denote 

/ -c 6 \ 
{) := I c —a , vu = V X u Wu. 
\-b a / 

The Poisson structure matrix is defined as the block diagonal 



J{z) 



\ 



Z12 



\ 



(4) 



Now becomes 
when we choose the Hamiltonian H 



z = J{z)VH{z) 



„ 1 y ^ Zij • Zab ~\~ 

ij {a,b)&NN{ij) i,j 



(5) 



where NN refers to "nearest neighbours" : 

NN{ij) = {zij-i, Zij^i, Zi_ij, Zi^ij}, 

and Hq represents the boundaries. For example, if we have zero boundaries {ziQ = 0, zqj = 
0, Zi^n+i = 0, Zn+i,j = 0), then 

Ho := 0, 

while if we have periodic boundary conditions, then 

Hq := ZQj ■ Zij + Zio ■ Zii- (6) 

We can easily extend this to a model covering both ferromagnet and antiferromagnet case. 



\ i,j NN i,j / 

where jx is the so called exchange integral [T|, assumed constant here, as in [2], and 

> for ferro 



(7) 



JK 



< for antiferro. 



Hence we have the Poisson system: 

z = J{z)VH{z), as in dH). 



(8) 



For an individual spin at the lattice point this becomes 

Zij = -jKZij X I ^ 2;afe - jKZij X Dz, 
\(a,b)eNN{ij) J 

= Zij X VH{z), 



(9) 



in both periodic and non-periodic cases. From now on we employ the notation 

NN{ij) {a,b)£NN{ij) 

Lemma 2.1. Any system of the form z = J{z)v{z) with ^ and v an arbitrary vector 
function, conserves the spin lengths in time: 



\zi,it)\ = \zij{o)\ yt,j,t. (10) 



Proof. 

d 



□ 

This gives us useful freedom in modelling. Next, the anisotropy term DS is approxi- 
mated by an average: 

Dzij D-{zi j^i + Zij+i + Zi_ij + Zi^ij). (11) 

this is sometimes referred to as the Roberts discretization. Now (jHl) becomes 

Zij = —jxZij X M(zjj_i -|- Zij+i + Zi-ij + ZiJ^ij) 

= Zi,xVH{z), ^ ^ 

where M = / + Z}/4isa diagonal matrix and H is modified according to (jll|) . 
Numerical method 

As we noted above, ((HI) is a Lie-Poisson system whose meaning we recall here: we can 
define 

{f,g}{z):=Vf{z)-{.J{z)Vgiz)), 
which fulfills the Jacobi identity 

{{f,9},h} + {{g,h}J} + {{hJ},g} = 0, 

hence {■, ■} is a Poisson bracket and J is a Poisson structure matrix. Since J is linear with 
respect to z, this Poisson structure can be derived from a Lie algebra structure, hence it 
is called a Lie-Poisson structure. 

For a detailed discussion on how to integrate this, see jTH] . To summarize that paper, 
the best way to integrate is to split the vector field in even-odd (or red-black) way: 



Zij = Vi + V2, 



(13) 



where 

{-jxZij X M EiviV(ii) even 

0, i + j odd, 

{0, i + j even 

X M 2;, « + j odd. 

Now, both of these flows can be exphcitly solved. For example Vi: for i + j odd Zij{t) = 
Zij{0). For i + j even, the sum over NN{ij) includes only pairs a,b with a + b odd, 
hence they are constants (during Vi). Likewise in V2 the sum is a constant. Denote the 
integrator of Vi by and that of V2 by $2,*. That is, 

<l>i,t = exp{tVi), $2,t = exp{tV2). 

The implemented integrator is a symmetric composition of these exact flows: 

:=$2,|0$i^iO$2 |. (14) 

This integrator 

• is time reversible 

• conserves spin lengths 

• in isotropic case {D = I) preserves energy 
since $i,t and $2,t do. See also Section ITTTl 



3 Dissipated version 

It is customary to add a dissipation term to ((H). In our case the corresponding dissipated 
version is derived from (fT^ and becomes 

Zij = Zij X VH{z) + azij X Zij x VH{z), (15) 

where a is a dissipation constant and the corresponding term is known as the Gilbert 
damping term. 

Clearly ()15|) can be written more compactly 

z = JVH + aJ^VH. (16) 

From lemma im it follows that 1^1 = 1 everywhere, i.e. the dissipation does not affect spin 
lengths. Let us first look at the Gilbert damping term more closely through the equation 

z = az X {z X B), (17) 

where 5; G M^, and a G M and B G are constants. Or, more compactly, 

z = aJ^B. 



This can be explicitly solved. Put 

V := z ■ B 



then (HH) is 



5 

w := z X B, 



V = a{—Ci + V 
w = avw, 
z = az X w. 



where Ci > constant, 

Ci = \z\^\B\^. 
We can solve for v (we have assumed \z(0)\ = 1): 



where 



S := exp{a\B\t), 
C2 :-- 



\B\ 


- vo 


\B\ 


+ Vo 



Note: if t — > oo, then 

a > => v(t) ^ —\B\ =^ z, B become antiparallel 
a < =^ v{t) \B\ =^ z,B become parallel. 

Substituting v we can solve for w, which is a scalar function times a constant vector: 
w{t) = f{t)w{0), 

= — Oast^oo,if«^0. 

Substituting w we can solve for z: 

z{t) = exp{gwo)z{0), 

, I IN sinf^fliyol) 

= cos (5( I Wo 1)^0 H . — . — X zo, 

I Wo I 

where 

t 

/C^ + 1 
j{r)dr = (arctanC — arctan(C£^)) , 

\B\C 



C := 





- Vo 




+ Vo 



Note that the exp above is a matrix exponential, while the sin and cos are the usual scalar 
functions. Here exp^gwo) is expanded as a Magnus series [TT]: the direction of w(t) is 
constant, hence qwq commutes with its integrals and Magnus series truncates after the 
first term. The evaluation of that term is by Rodriguez' formula, hence (jSBI)- 

Evaluating g numerically was a problem because eventually v approaches ±|-B| (physi- 
cally this means z becomes (anti-)parallel to B) so the C in g becomes zero, g itself is not 
singular, however this presentation is difficult to evaluate. We used the following Taylor 
expansions in the implementation: if — fo| < 0.0001, 

g = -l + £ + C'(^-^ + £- y-"^ + 0{C'), 

and if \\B\ + vo\ < 0.0001, 

g = l-S~' + C-'(l-S-' + ls~A + 0{C-% 



A Lyapunov function 

Note that 

—B i) 

\w\^ = {z X B) ■ {z X B) = -{B X z) ■ {z X B) = -B ■ z X {z X B) = z = , 

a a 

hence 

= -a\z X 5p < 0, if a > 0. 

dt ' I - ' - 

So f is a Lyapunov function, when a is positive. If H{z) := Cv = Cz ■ B, C constant 
scalar, then 

= -aC\zx B\\ (26) 

that is, if is a lyapunov function iff sgn(Q;C) = 1. Later sgn(C) chooses between ferro- 
magnet and antiferromagnet. 



Several spins 

Now we continue from (fT3j) . which can be written 

Zij = Zij X {-Jk)M ^ z + azij X Zij x {-Jk)M ^ z. (27) 

NN{ij) NN{ij) 

Recall from the previous discussion that a < implies Zij tends to become parallel to 

i-jK)M J2 ^■ 

NN{ij) 

This means, see fEHjl . that if jx > 0, then sgn^ajx) = —1 and energy H is decreasing. In 
other words, for a ferromagnet negative a means energy damping. 

To summarize, a ferromagnetic or antiferromagnetic spin system subject to Gilbert 
damping will uniformly dissipate energy for appropriate choice of the sign of the damping 
coefficient. Moreover, a Gilbert-damped system is spin-length conserving. 



Numerical method 



To integrate, we split the vector field in even-odd way as in the conservative case (Section 

El 

% = Vi + V2 + V^ + \/4, (28) 



where 



-3KZ^j X M Y.NN{ij) ^ + j even 

i + j odd, 

0, i + j even 

-jxZij X Mj^NNiij) ^ + i odd. 



^ ^ 1 -Jxazij X Zij X MY,NN{ij) ^ + i even 

i + j odd. 



V^4 



, i + j even 

-jxazij X X M X;jv7v(ij) 2;, i + j odd. 



Now, all these fiows can be explicitly solved. For example Vi. for i + j odd %(t) = %(0). 
For i + j even, the sum over NN{ij) includes only pairs a, b with a + b odd, hence they 
are constants (during Vi). Likewise in V2, V3, and V4 the sums include only constants. 
Hence in Vi and V2 we solve 

Zij = Zij X B, B constant, (29) 

and in V3 and V4 we solve 

Zij = azij X Zij X B, B constant, (30) 

which are solved above. Note that (j^^ . ()30|) have different 5's. The implemented inte- 
grator is a symmetric composition of these exact flows: 

where = exp(t V^) are the exact flows. 

An important feature of our method is that it dissipates energy when the flow p5|) does. 
This can be seen in the following way: we solve the flows Vi, . . . , V4 exactly, hence every 
step in the composition pip follows the energy of the associated vector fleld exactly. Now 
$At is a second order method and it follows the energy evolution with accuracy O(At^). 
With small enough time step At the error is negligible and our method dissipates the 
energy. 



4 Thermostatted version 

The motivation behind using thermostats is keeping the system around some constant 
average temperature. This is a reasonable assumption for example in systems with heat 
baths. That is, we allow the energy to fluctuate. But, at the same time, we want to keep 
the spin lengths constant. This will be carried out by modifying the dissipation term 
introduced in previous sections. 



We use the thermostatting term suggested in 2j: choose a parameter T (temperature) 
and compare the system's energy to it, allow the damping coefficient a to vary with time: 
a = a{^) and 

« = - (}^)' E - ^^^-«.) ■ ^ % ^ ^) ' (32) 

where T is temperature, and k is Boltzmann's constant which we hereafter take to be 1. 
M is number of degrees of freedom, that is A/" = 'ir? since we have an n x n square lattice, 
K is coupling strength and typically ~ \fM . Here we take KjM := 1/n. X is equal to: 

X = V^^M := -j^M ^ z. (33) 

The thermostatting variable a has been given the nickname "global demon" , so called 
due to its non-local (hence non-physical) character: it affects all spins simultaneously. 
So our thermostatted system is 

Zij = Zij X X + azij X Zij X X (34) 
« = -(-^)'$^(X-TV.J-(z., xz,, xX). (35) 

It is possible to show that the ferromagnetic system thermostatted using (jMjl . ljH^ 
samples from the canonical ensemble. This system also conserves spin length. Finally, 
one easily demonstrates that these equations are invariant under the simultaneous time- 
coordinate transformation t \—>- —t, z t-^ —z^ a ^— > —a, i.e. the equations are time- 
reversible. 



Numerical method 

To integrate, we split the vector field in even-odd way as above, with the a term: 



Zij 



a 



V^l + V^2 + V^3 + V^4 + V^5, (36) 



where V^i, V2, V3, V4 as in Section El (with a = =constant) and 

Zj,- = constant, 
^ { (37) 
« = E., ((% -Bf-B-B- 2Tz,, . B) . 



Here we have simplified: 



B := X = — jx M z = indep. of 

NN{ij) 



Zij, 



V,- {z X z X B) = V,- {{z ■ B)z -B) = 2z- B, 

B ■ {zx zx B) = {z- Bf - B ■ B. 

In V5 all terms are constant so the equation with V5 is trivially solved. But note that 
the update step of a is Oijn?). 

:= $14 o $2 * o $3 t o $4 . o $5,4 o $4 * o $3 * o $2 . o (38) 



where $j^f = exp(t V^) are the exact flows. 

An important feature of this discretization is that it is time-reversible with respect to 
the mapping z i— >■ —2;, a ^— —a, t ^ —t. This can be seen by recalling from Section 
II .11 that if f{Ru) = —Rf{u) for some linear involution R, then the Strang splitting gives 
a time-reversible map. Here u := {z,a) and Ru := {—z,—a). Applying the rule four 
times in a row: first to $4 and $5 in the roles of $t,/i and $r,/2 of Section secondly 
to $3 and $4 o $5 in a similar way, next to $2 and $3 o $4 o $5 and finally to $1 and 
$2 o $3 o $4 o $5, we get the claim. 

5 Numerical results 

In all our simulations we used n = 50, that is, a 50 x 50 lattice. We used ferromagnets with 
anisotropy: D =diag(l, 1, A). This is known as "easy plane" or "easy axis" anisotropy, 
corresponding to A < 1 or A > 1, respectively. 

5.1 Dissipated system 

Example 1. In Figure^ we see evidence of the dissipation of energy. On the top part is 
energy in semilog scale, on the bottom part is the maximum norm of the discrete Laplacian 
during each time step. Here we used periodic boundary conditions, A = 1.1, a = —0.5, 
and timestep At = 0.1. The initial configuration was random (top left of Figure |21). 

The evolution of the discrete Laplacian is understood as the system settling down to 
some formation, and this can be seen in Figure El which includes snapshots of the same 
simulation. The snapshots describe the components of the spins. The order of the 
pictures is top row first, from left to right. The darker a point is, the lower component 
it has. Black represents spin down, white spin up. 

The result is typical: the dissipated system converges to two bands of up and down 
spins. We also tested zero boundary conditions, then the dissipated system typically 
converged to a single band. 

5.2 Thermostatted system 

Example 2. In Figures El and 0] is a thermostatted system with periodic boundaries, 
A = 0.9, T = 0.04 and timestep At = 0.01. The initial condition is random. In the top 
part of Figure El is the thermostatting variable a and in the middle part is the energy, 
and in the bottom part the maximum norm of the discrete Laplacian. After an initial 
phase both a and energy settle to an aperiodic oscillatory motion, a between —10 and -|-10, 
energy between —4800 and —4300. We plotted only 2000 steps but the behavior continued 
similarly for at least 25000 steps. In Figure 01 we can see slowly creeping boundaries; the 
reader is asked to compare the white areas. At t = 0.26 it suddenly looks chaotic, then 
renders back to the creeping boundaries. This can be seen as a kind of stability of the 
creeping boundaries. After 25000 steps there still is slow motion, the system does not 
converge to any particular formation. 

Example 3. Another example, in Figures OHl is a thermostatted system we call "wan- 
dering vortices" . This beautiful system has random initial conditions, periodic boundaries, 
and parameters A = 0.9, T = 0.05 and timestep At = 0.05. Figures El and El represent the 
evolution of a, energy, and averages of the energy over different time windows. Interest- 
ingly, the behavior of a is much wilder than in the Example 2. The snapshots in Figures 
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Figure 1: Energy of dissipated system. 



I7p8l show the 2;— components of the lattice. From random state, the system very quickly 
forms vortices on smooth surrounding, which move around for a short time, then look 
random again, then vortices again. Sometimes the vortices died out completely leaving 
us just with smooth surface, then reappeared again. 

Comparison to RK4 with projection 

For comparison we implemented the classical Runge-Kutta 4*^ order method (RK4) with 
projection: after every step we normalize 

where Zij^Rx denotes the result of RK4 step. At very small timesteps for which the RK4 
method could successfully integrate the problem, it was slightly faster than splitting, 
but the RK4 method became rapidly unstable as the stepsize and/or anisotropy were 
increased. The sphtting method was able to handle large anisotropics (A = 3) and step 
sizes (At = 0.3). However, we did not seek the limits of our splitting method. The values 
A = 3 and At = 0.3 indicate the superior stability well enough at this stage. 

We tested this projected RK4 on Examples 2 and 3. We kept the other parameters 
intact but changed the step size. As a sign of failure, we stopped the computation when 




Figure 2: Snapshots of 2;-components (black is down spin, white is up spin) of the dissi- 
pated system, example 1. 
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Figure 3: Energy, alpha, and norm of Laplacian of thermostatted system. 



the code started to produce infinities. In Example 2 the maximal timestep was 0.01, while 
in Example 3 At^aa; = 0.015. The results are summarized in tabled 



6 Discussion 



In this paper we have developed and tested a geometric integrator for a semi-discretized 
Landau-Lifshitz-Gilbert (LLG) equation which includes a nonlinear dissipative term, as 
well as for a more complicated thermostatted model, following the approach of Bulgac 
and Kusnetsov. 

The integrator for the dissipated system is shown to have a dissipative property. How- 
ever, it is difficult to compare since we do not know the exact continuous solution. LLG is 
currently a very active topic of research, see more details in the introduction. However, it 
seems that so far there has not been developed a geometric integrator for a thermostatted 
system. 

Trying to simulate the thermostatted systems with projected RK4 revealed both the 
features of a stiff ODE and features of a conservative system. The combination is ex- 
tremely difficult for standard form numerical methods. The key feature of our splitting 
method is that it is constructed from composition of building blocks that simulate each 
of the two components of the system correctly. 




Figure 4: Snapshots of ^-components (black is down spin, white is up spin) of thermostat- 
ted system, Example 2. 
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Figure 5: Energy and alpha of "wandering vortices", Example 3. 



Simulation with our new thermostatted method has revealed interesting phenomena: 
slowly creeping boundaries, or wandering vortices, both of which appear from random 
initial conditions. Our informal term "wandering vortices" in Example 3 refers not to 
certain particular vortices that survive throughout the whole simulation, but to a situation 
where we have two or more vortices which wander for a while, then violently crash and 
form new vortices. Intermediate states include "quasi-chaotic" state, an informal term 
by which we mean a state that suddenly appears and looks random but is not, since it 
renders immediately back to (almost) the same smooth motion. 

The RK4 method is less stable. The stepsize restriction is an order of magnitude 
smaller compared to our splitting method. This is evidence of stiffness in the ODEs, and 
a better choice might seem to be a stiff solver on this account, but if one uses a stiff 
solver the result would be poor resolution of the conservative evolution which is also an 
important component of the dynamics of the system. The best compromise is therefore 
a composition scheme, such as that outlined here, which separately and appropriately 
resolves each term of the system. 

We anticipate that this work will stimulate further research in the development of 
thermostatted numerical methods for systems with complicated nonlinear structure." 
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Figure 6: Averages of energy of "wandering vortices", Example 3. 
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Example 3 



At 
0.020 
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Table 1: Failure of projected RK4 
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Appendix: Lie-Poisson Canonical Sampling Technique 

Consider first the Poisson rigid body system consisting of a Hamiltonian H{z) together 
with structure matrix J{z) admitting Casimir \z\. Based partly on [12. we construct an 
augmented Hamiltonian with additional variables cr, tt„, 6, hq: 

H = H{Q{e)z) + 7r2/2/i + A;T In a + G{d, tt^, tt,). 

where Q is an orthogonal matrix depending on parameter(s) 9. The Lie-Poisson structure 
is just the rigid body Poisson structure augmented by the canonical structure for the 
augmenting variables. Under assumption of ergodicity, and some very mild technical 
conditions similar to those obtained in ^2], this Hamiltonian can be shown to provide 
canonical sampling from microcanonical trajectories i.e. 

f{Q{e)z)6[H - E]dadedn,d7re = f{z) exp{-^H{z)), 

with preservation of the Casimir due to Q being orthogonal. In order to obtain ergodicity, 
the "bath Hamiltonian" G should be sufficiently complicated. 

In the case of a spin system, H{zi, Z2, . . . , z^), we may introduce a separate unit 
3-vector 6i for each spin vector. Then the Hamiltonian 

H = H{Q{9i)z^, Qie2)z2, . . . , QieN)zN) + Trl/2fi + kT\na + 0(6, ne, vr,), 

will enable canonical sampling. For example Q{u) can be a Householder transformation, 




Q(u) = /-2 



T 

UU 



U 

and the bath Hamiltonian G can describe a coupled system of spherical pendula involving 
in some nontrivial way the parameter tt^. 



